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Abstract 

We apply the adiabatic time-dependent density functional theory to magnetic circular dichroism 
(MCD) spectra using the real-space, real-time computational method. The standard formulas 
for the MCD response and its A and B terms are derived from the observables in the time- 
dependent wave function. We find real time method is well suited for calculating the overall 
spectrum, particularly at higher excitation energies where individual excited states are numerous 
and overlapping. The MCD sum rules are derived and intepreted in the real-time formalism; we 
find that they are very useful for normalization purposes and assessing the accuracy of the theory. 
The method is applied to MCD spectrum of Cgo using the adiabatic energy functional from the 
local density approximation. The theory correctly predicts the signs of the A and B terms for the 
lowest allowed excitations. However, the magnitudes of the terms only show qualitative agreement 
with experiment. 
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I. INTRODUCTION 



Magnetic circular dichroism (MCD) is an important spectroscopic observable useful for 

n n 

characterizing the electronic structure of molecules [1[ and condensed matter systems 
The theory of the MCD response is usually given as a third-order perturbation in a basis 
that diagonalizes the zero-field Hamiltonian. This formulation, called the sum-over-states 
method, requires a considerable computational effort to construct the states and perform 
the summations. There have been recent attempts to simplify the calculation by using 
eigenstates of the Hamiltonian in the presence of the magnetic field , but one still 
requires a sum over transition densities. We propose here a completely different formalism 
based on the solution of time-dependent equations of motion and present a formalism for 
A and B terms of the MCD. We find that the formalism is a practical one when applied in the 
framework of time-dependent density functional theory (TDDFT). In fact, the TDDFT has 
already been used successfully to calculate MCD in molecules . A separate problem in 
the theory of MCD is the choice of a basis set to construct the electron orbital wave functions. 
The MCD puts higher demands on the orbital representation to satisfy completeness and 
gauge invariance. In our treatment, we represent the orbital wave functions on a spatial mesh 
rather than with atom-centered basis set. The calculated matrix elements are automatically 
gauge-invariant and one also avoids the inconsistencies that cause sum rules to be violated. 

uhe real-time TDDFT has been applied to many observables related to 
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id]. Specific applications include the molecular absorption spectrum 



hyperpolarizabilities 12Nl4|. the dielectric function [15], and chiral 



electron dynamics 
in the continuum 

dichroism [l^. The real-time method has also been applied to phenonema associated with 
high fields. In the presence of high fields, there is hardly any alternative theory available, 
at least at the ab initio DFT level. Applications include nonlinear electron dynamics in 



metallic clusters 



17l |. high harmonic generation [l8l. Coulomb explosion 
breakdown 2l|], and coherent phonon generation 221]. 
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The organization of this article is as follows. In Section II we define the time-dependent 
quantities that are computed and derive the formulas for extracting the observables related 
to MCD. We also review the sum rules satisfied by the MCD response in that Section. In 
Section III we provide some of the numerical details in carrying out the MCD calculations. 
Following that, in Section IV, we apply the theory to the Cqo molecule. Due to its high 
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symmetry, the Ceo molecule can exhibit both A and B terms of MCD. The lowest electronic 
excitations of this molecule are the tt — tt* character; for the measured transitions we find 
the correct signs for the calculated A and B terms. The theory is in quahtative agreement 
also with magnitude of the B term of the lowest transition. However, the present theory 
does not reproduce well the other transitions and the magnitude of the A term. 

II. THEORY 

A. Definitions for MCD response function 

We consider a molecule under a static magnetic field S in 7 direction. The electronic 
Hamiltonian is written as 

H = Ho + iibBL^, (1) 

where Hq is the Hamiltonian in the absence of the magnetic field, is the angular momen- 
tum operator, and /ib — e/2mc is the Bohr magneton. We take a convention of e > and 
Ti — l. We denote the ground and excited states under the static magnetic field B as 

= EA. (2) 

We denote the dipole operator as — — eX^f, where rj are electron coordinates. We 
define the circularly polarized dipole operators with the normalization convention = 
{Ha ± ^A*/3)/'\/2, where (q;^7) is a cyclic permutation of (xyz). 

In MCD, the basic object of study is the difference in dipole absorption strength functions 
for hght of opposite circular polarization in a weak magnetic field. The MCD response may 
be defined by the strength function 

RmCb{E) = ;^EE'^(^-^n0){|(<f0|/xL^^|<f„)r- |($0|/xt^|$n)|n, (3) 

where £'„o is the excitation energy of state n, E^q — E^ — Eq. The beam direction coincides 
with the magnetic field direction along an axis labeled by 7. It is convenient to express the 
strength function in a Cartesian basis using the antisymmetric tensor Ca/s-y, 

3//Bi^ a/37 

where the magnetic field is applied to 7 direction. 
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To assess the quality of the theory, it is also useful to calculate the ordinary dipole 
response. We define the dipole response Rd{E) 

Rn{E) = Wj2S{E-E,,omo\f^a\^n)\', (5) 

^ n a 

This is related to the oscillator strength distribution by 

df 2mEno 



dE e2 



Rd{E). (6) 



Below the ionization threshold, electronic states are discrete. In this energy region, the 
MCD strength function is often written as 



) , (7) 



where Tn{E) is the spectral shape of the n-th state normalized as / Fn{E)dE = 1. The 
zero-field excitation energy is expressed as e\^~^\ There appear both A and B terms for 
molecules with degeneracy in either ground or excited states, while only B term appears for 
molecules without degeneracy in any states. Integrating the MCD response function over 
an excitation energy around E^~^\ we have 

An = /^"°^' dE{E - eS=''^)Rmcd{E), (8) 

and 

rEno+e 

Bn = / dERucBiE), (9) 

where e is a small energy interval. 

Employing the perturbation theory, these terms may be expressed in terms of the energy 
and the wave functions in the absence of the magnetic field. For A, we have 

An = -\Yeo,)3jY.{i^l)nt,nt ' (L^)o.,0 J Ini($Os | $nt) ($nt I)" I *0s) , (10) 

a/37 

where s and t distinguishes degenerate states of ground and excited states, respectively. The 
basis which diagonalize is assumed. For B, we have 

2 f 1 

Bn = -Im^e„/37 Hi T^(*m|^7l^o)($o|/Wa|4'n)(^'n|/^/3|$m) 
•J a/37 ^^mO 

+ -^($n|^7l*rn)(<I'o|Aia|$n)($m|/^/3|$o) j ■ (H) 
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Similarly, we define the ordinary dipole strength as 

Vn = j''^'^^ dERoiE) = |($o|/i„|$„)|^ (12) 

It is related to the oscillator strength /„ as /„ = 2mEno'Dn/ . 

Finally, with our definition of the MCD response, the B coeffient is related to the ratio 
of the MCD extinction coefficient Ae to the ordinary extinction coefficient e by the formula 



e'fn 



(13) 



B. Real time formulation 



The response associated with any pair of operators, Oi and O2, may be calculated in 
the real time domain starting from the ground state wave function $o- One first applies 
an impulsive perturbation XOi to obtain an initial wave function \E'(t = 0+). This is then 
evolved in time with the exponentiated Hamiltonian operator, 

^(t) = e-'^V^^^$o. (14) 

The real-time response S2i(t) is given by the expectation value of the operator O2 in that 

state, 

S2l{t) = {^{t)\02\^{t)). (15) 

The linear response is evaluated by treating A as a small parameter and taking the deriva- 
tive dS2i{t)/d\ at A = 0. Depending on the operators Oi and 02-, the strength function 
R21 is obtained as the imaginary or real part of the Fourier transform of linear response 
dS2i{t) / d\\\=Q on the time interval [0, +00]. 

This general formulation of the linear response applies to the MCD strength function Eq. 

(jl]) taking the operators to be Oi = and O2 = Ha- The wave function is set up at t = 

7) 

S^:}{t) = {^>,f,{t)\i,^\^>up{t)), (16) 



as \l/A;/3(t = 0+) = e*'^^''$o and the real-time response S^}{t) is given by 



where (7) in S^^/^ is included to remember that a static magnetic field is applied to 7 direction 
throughout the time evolution. Expanding the perturbing operator e^^'^i^ in powers of k, we 
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have 

= zA;$]($o|/Xa|$n)($n|/X/3|$o)e-^^"«*-^fc5](<l'o|/i/3|$„)($n|/ia|$o)e'''"''*. (17) 

n n 

In the last formula, we have expressed the time-evolution operator in terms of the energy 
eigenstates of the system. We next separate out the time-even and time-odd parts of the 
response, writing it as 

S^J^it) = -2A;^Im {(<l>o|/x„|$„)($„|/i/3|$o)}cosKot 

n 

+2kJ2^e {($o|Aia|$n)($n|/X/3|*o)}sinE„ot. (18) 

n 

The combination S^J^ - S^J^ = -AkJ^rJ^^ {{%\l^a\^n) {^n\lJ'fi\^o)} cos Enot 
isolates the first term with the cosine dependence on time. We obtain an expression propor- 
tional to the MCD response by taking its cosine Fourier transform, 

- A rdtcosEt{S^J^ - S'^;n = J^Im {{<l>o\fi^\<l>n){^n\fip\%)} SiE - E^o)- (19) 
Zttk Jo n 

The MCD strength function with the proper prefactor is given by the integral over the 

real-time response as 

Rmcb{E) = E r dtcosEt {S^^} - ^g} . (20) 

For the molecules we treat here, we can choose coordinate systems such that the off-diagonal 
response is antisymmetric, S^^{t) = —Sj^^{t). Then the second term in Eq. flTK]) is identi- 
cally zero and Eq. fl20l) reduces to 

Rmcd{E) = f dtcosEt {Si'jit) + 4:)(t) + SiyJit)} . (21) 

This is our main result that will be applied to calculate the MCD. 

For most if not all MCD spectra, the sign of Rmcd{E) on the infrared side of the lowest 
optical excitation is negative. We shall call this the "normal" sign. 

The ordinary dipole response Rd{E) may also be computed in the formalism as Fourier 
sine transform, 

Rb{E) = ^Y.I dtsmEtS^^2it)- (22) 

Note that this can easily be evaluated at the same time as Rmcd{E) because the same 
time- dependent wave function is used in both. 
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C. Sum rules 



The real-time formalism is very convenient for evaluating and verifying energy-weighted 
sum rules. In particular, the MCD response satisfies a quadratic sum rule that depends 



only on the magnetic field strength and the number of electrons 23|]. The connection to the 
til 
t, 



time-dependent response may be easily derived by expanding S^/^ as a power series in time 



S;^^{t)c^so + S2t' + ---. (23) 

Only even powers of t are present in the series expansion, due to the suppression of the 
second term in Eq. f|T8l) . The coeffients sq and S2 can be readily expressed as commutators 
of the /i operators and the Hamiltonian and evaluated analytically. One finds 

So = -2/cIm^{($o|/ia|$n)($n|/i/3|$o)} = ik{<l>o\[fla, = 0, (24) 

n 

and 

S2 = fcIm5:{($o|/i.|$„)($n|/i/3|$o)}i?^0 = -y(*0|[/ia, [H, [H,^!^^) = ±^^^e 

' (25) 

where A^e is the number of electrons and the sign reflects the order of (a/37). These formulas 
can be expressed as integrals over the MCD strength function 

dERucDiE) = (26) 



and 

POO 26^ 

I2 = / dEE^RucDiE) = —N,. (27) 
JO m 

The I2 sum rule has a simple physical interpretation. Consider the real-time response 
associated with S^^J. The impulsive exciting potential eky 5{t) gives the electrons an average 
momentum equal to —eky, the integral of the force over time. The corresponding velocity, 
V = —eky/m, is subject to a magnetic force —ev x B/c which is in the x direction for our 
geometry. The corresponding acceleration is a = e^Bkx / m?c. Thus the x-component of the 
dipole moment increases quadratically with time according to the acceleration formula 

(— (*» = -,at' = (28) 

in agreement with Eq. (!23|) and (!25|l . 



In the results of the calculations given below, we will also show the performance of the 
theory with respect to the /-sum rule. In terms of the response Rd{E), we define the 
quantity 

/e = 2m r dE'E'RD{E'). (29) 

The sum rule is given by 

/+00 = (30) 

where is the number of electrons. The associated short-time behaviour of the real-time 
response is linear in t and given by 

{^kait)\li^\^kait)) ^ N,t. (31) 

m 

D. Time-dependent density functional theory 

The basic variables in Kohn-Sham density functional theory are the orbitals (piir), which 
are varied to minimize an energy expression E^s that contains the quantum kinetic operator 
and terms depending on the density n(r*) = J2i The formal variation of the E^s 

energy expression with respect to an orbital gives the Kohn-Sham operator H^s- In the 
time-dependent extension of DFT with the adiabatic approximation, the operator Hks also 
used in the equation of motion for the orbitals. To linear order in the magnetic field strength 
B, the time- dependent orbitals satisfy the equations 

HKs + I^BL-B\^i{t)=t—^i{t). (32) 

These equations are solved for ipi{t) with initial condition tpi{t = 0^) = e^^^f' (j)i{r) . where (pi 
are the ground-state Kohn-Sham orbitals. The time-dependent response is calculated as 

S^:^{t) = Y.{Mt)M^{t)). (33) 

i 

For the calculations described below, we treat only the valence electrons dynamically and ig- 
nore any spin dependence. The effects of core electrons are treated by using Troullier-Martins 
pseudopotentials for the electron- ion interaction [2J]. We use the usual LDA functional [25 1 



for the exchange-correlation interaction as in previous work 



26|. 
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III. COMPUTATIONAL ASPECTS 

The implementation of our real-time TDDFT is described in detail in Ref. ^|. An im- 
portant difference from other TDDFT codes is the orbital representation in 3-dimensional 
Cartesian mesh. This has the computational benefit that the Kohn-Sham operator is given 
by a sparse matrix. The representation is also convenient for treating extended wave func- 
tions such as Rydberg states or continuum states. It has the disadvantage, however, that it 
does not permit all-electron calculations with practical mesh sizes. For checking our code, 
we found it helpful to calculate the observables with the Troullier-Martins pseudopotential 
replaced by an anisotropic harmonic oscillator potential. All the observables in this model 
have analytic expressions that can be compared with the calculated numerical quantitites. 
See the Appendix for details. 

In our implementation of the mesh representation, the momentum operator p is computed 
by the 8-point difference function, which is consistent with our treatment of the kinetic 
operator /2m as a sum of 9-point difference functions along the three Cartesian axes. 
The main numerical parameters in the calculation are the mesh spacing Ax, the size of the 
spatial domain on which the orbital wave functions are defined, the time step At, and the 
total integration time T. We find that adequate precision for our purposes is obtained with 
parameter values Ax = 0.5 and At = 0.03 in atomic units. The spatial dimensions needed 
for the orbital wave functions depend on the desired accuracy in the continuum region. The 
continuum strength functions are smooth only if the spatial domain is large and absorbing 
boundary conditions are applied at the edges. Typically, we take a cubical box of 160'^ mesh 
points for the calculations. For small molecules, a much smaller domain is adequate if the 
details of the response in the continuum are not needed. 

Although the TDDFT is fundamentally nonperturbative, the quantities we calculate are 
in fact the perturbative limits with respect to the strengths of the applied magnetic and 
electric fields. We thus choose strengths that are small enough for the linear response 
formula to apply, but large enough to avoid numerical roundoff errors. For the perturbing 
electric field, we take k = 0.001. For the magnetic field, the calculations reported below were 
carried out with a magnetic field given by ^bB = 0.0005 au. The intensity of this magnetic 
field is 0.137 au. For comparison, a field strength of one Telsa has the value 5.81 x 10~^ in 
atomic units. 
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The integration time T required to calculate the response depends on the desired energy 
resolution. We multiply the integrands in the Fourier transforms Eq. f l20|) and Eq. f l22|) by 
the filter function 1 - 3(t/T)2 + 2{t/Tf to smooth out spurious oscillations from the upper 
time cutoff. The resulting peaks associated with sharp states have a width F (full width at 
half maximum) given approximately by F 6/T. Most of our results were calculated by 
integrating Nt = 60000 time steps, giving F ^ 6/{NtAt) ~ 0.0033 au = 0.1 eV. 



IV. APPLICATION TO Ceo 

The Ceo molecule offers a good test of the methodology to demonstrate the feasibilty of 
using the real-time method as applied to fairly large molecules. Due to the high symmetry of 
Ceo, all optically allowed transitions are three-fold degenerate and there will be both A and 
B terms in the MCD spectrum. There are 5-6 excitations in the calculated spectrum up to 
6 eV, all of which are vr — vr* character. It has been found that the experimental oscillator 
strength 27, [28I accords well with the theory 26 1 based on the ab initio adiabatic local 



density approximation. 



26| for the oscillator 



Our calculation here is very similar to that carried out in Ref. 
strength function. The integration time in the present calculation is somewhat longer, 60,000 
time steps with At = 0.03 au. Before examining the MCD response, we recall the results for 
the ordinary dipole response, as calculated in the real-time method. Fig. [T] shows the S^zit) 
real-time response over the interval [0,T] = [0,25] fs with the left-hand panel showing an 
expanded view of the first 0.275 fs time interval. The short-time behavior expected from Eq. 
(!3T|) is shown by the straight dotted line in the left-hand panel. One may see that the initial 
response does indeed follow Eq. fl3Tl) very well. After the initial rise in the first 0.1 fs the 
dipole moment oscillates with a period of order of one fs correspond to the strong transitions 
in the energy interval 7-15 eV. Note that the oscillation is essentially undamped. This is 
a consequence of the sharpness of the bound excitations that would produce a 5-function 
response if the Fourier transform could be done exactly. The numerical Fourier transform 
was carried out to final time t = 1800 au = 43.5 fs with the results for the low-frequency part 
of Rd{E) shown in the left-hand panel of Fig. [2l There are four transitions in the spectral 
region 0-6 eV, at excitation energies of 3.5,4.3,5.3 and 5.9 eV. The numerical FWHM widths 
are about 0.1 eV, as expected from the integration time. The important information besides 
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FIG. 1: Real-time dipole response Szzif) = {'^ kz{t)\fji'z\^ kz{t)) for Ceo- The left-hand panel shows 
the time interval t = — 0.275 fs with a linear time scale. The sloping line shows the expected short- 
time behavior according to Eq. (j3ip . The right-hand panel shows the time interval t = 0.25 — 25.0 
fs on a logarithmic time scale. 
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FIG. 2: Dipole response for Cqo- The differential oscillator strength df/dE (Eq. ([6])) up to 6 
eV is shown in the left-hand panel. The right-hand panel shows the integrated oscillator strength 
function Je, Eq. ([29]) . 

the transition energy is total strength in the individual peaks. This can be extracted from 
the graph of the integrated strength defined in Eq. ( l29ll . The jumps at low energies give 
the / strengths of the discrete transitions. The total integrated strength is / = 233, rather 
close to the sum rule number / = 240 for the N^. = 240 valence electrons treated dynamically 
in the TDDFT. We note that the sum rule is not expected to be satisfied exactly for our 
energy functional, because of nonlocality in the Troullier-Martins pseudopotential. 
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FIG. 3: MOD real-time response Sxy in Cqo- The solid line in left-hand panel shows the evolution 
for short times, < t < 0.3 fs. The dotted curve is the expected dependence from Eq. (j28p . The 
long-dashed curve shows the response in the time range < t < 0.07 fs with the nonlocality in the 
pseudopotential turned off. The right-hand panel shows the response in the longer time interval 
0.25 < t < 40 fs on a logarithmic time scale and a magnified ordinate scale. 

We now take up the MCD response. The left-hand panel of Fig. 3 shows the calculated 
MCD real-time response S^^{t) over the time interval < t < 0.3 fs. The dashed curve in 
the left-hand panel shows the predicted short-time dependence according to Eq. (!28|) . The 
computed time dependence starts out quadratic as expected, but the coefficient of is lower 
by 40% than expected from Eq. ( |28l) . To confirm that the nonlocality of the pseudopotential 
is responsible for the disagreement, we have recomputed the response for short times with 
nonlocality of the pseudopotentials turned off, shown as the long-dashed line in the Figure. 
This agrees closely with the expected short-time behavior. We do not have any explanation 
why the sum rule violation is much stronger for the MCD strength than for the ordinary 
dipole strength. 

The MCD response going to long times is shown on the right-hand panel of Fig. |3l It 
is interesting to note that the amplitude of oscillation increases with time. This behavior 
is in contrast to the ordinary dipole response, which has a maximum excursion in the first 
oscillation. The reason for the increase in amplitude is the presence of the A terms which 
give a real-time response that cancels at short times and only becomes visible at later 
times. Taking the Fourier cosine transform of the real-time response using Eq. ( l20i) . we 
find the MCD spectrum shown in Fig. HI left-hand panel. The ^-type character of the 
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7T—7T* transitions is clearly seen in the shape of curves, each with a strong alternation of sign 
over the width of the peaks in the dipole response function. (Again, there is no physical 
significance to the calculated widths since they depend on the integration limit in the Fourier 
transform.) It is interesting to see that the sign of A coefficients can vary from state to state. 
The excitation at 5.9 eV has the normal sign, namely negative on the low-frequency side, 
but the three lower excitations have the opposite sign. The four transitions in the figure 
also have a significant B-tjpe MCD response, visible by unequal positive- and negative-going 
peaks on the two sides of the transition. The i3-type response may be seen more clearly 
in the graph of the integrated MCD response, dE'RucDiE'), shown on the right-hand 
panel of Fig. m The Bn coefficients can be read off from the step increases going across each 
transition, cf. Eq. (9). The values are reported in Table I, divided by the theoretical dipoles 
strengths Vn (Eq. (12)). This is to facilitate comparison to the experimental values 29|, 
which are given in this form. We see that the signs of Bn for the lowest two states agree. 
This is far from trivial. Also, the calculated magnitude of the lower one is within a factor 
of 2 of experiment. This is poorer agreement than is typical for the calculation of oscillator 
strength /„ in TDDFT, but perhaps this should not be unexpected due to the difficulties 
uncovered by the unexpected short-time behavior. Also, we know that there is considerable 
screening of the valence electron transition moments, amplifying the relative errors of the 
screened observables. The Bn of the second state has a much larger discrepency. Until that 
is understood, one cannot use the TDDFT as a predictive tool for large molecules. 

For a overall view of the MCD response, the right-hand panel of Fig. H] shows the inte- 
grated MCD response up to 50 eV. The integrated response is predominantly negative, as to 
be expected with the negative-going initial evolution. One sees that the total goes to zero 
at the upper energy, showing that the Iq sum rule (Eq. (26)) is nearly satisfied. Finally, the 
E"^ sum rule, Eq. ( 127|) . has a value I2 = 258, almost a factor of two smaller than the nominal 
value of 2Ne = 480. We have already seen this effect of the nonlocality in the short-time 
response. 

We next turn the ^-type response, arising from the energy splitting between members of 
^Ti„ multiplets as in the Zeeman splitting. A convenient way to express the splitting is as 
the effective (^-factor for the transition |l][Eq. (52)], 
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FIG. 4: MOD response RyiCDiE) in Cgo- Left-hand panel shows the strength function Eq. Q. 
The corresponding integrated strength function is shown in the right-hand panel. 

TABLE L MCD response of the lowest four ^Ti^j states in Ceo- The experimental data is from Ref. 
29( 1 . Our calculations are given in the columns labelled TDDFT. The effective orbital 5- factor is 
defined in Eq. ([31 



Energy (eV) 
Exp. TDDFT 


Exp. TDDFT 


9 

Exp. TDDFT Ref. [30] 


3.8 3.5 

4.9 4.3 
6.0 5.3 

5.9 


100 64 
-700 -146 
66 
-120 


-0.3 ±0.05 -0.97 -1.0 
-0.55 ±0.15 -0.58 -0.75 

-0.20 ±0.12 

±0.35 



This is related to An by 



An 

2pr 



(35) 



We extract the An coefficients by Eq. ([8]) from Rmcd{E). We may also extract the energy 
shift AE from the zero-field value using the formula 

, Rmcd (E) - Je"°^' dE -Rmcd {E) 



AE 



(36) 



4i?D(^no) 

The extracted (yf-factors are shown in Table I along with the measured values j29| and 
results of a model calculation jsol . As with the Bn values, we see agreement on sign for the 
two measured transitions. However, only the upper transition has a magnitude consistent 
with experiment. 
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V. CONCLUDING REMARKS 



We have shown that from a computational point of view, the real-time method is a practi- 
cal approach to calculate the MCD response in TDDFT. In particular, the entire response in 
the energy region of valence-electron excitations is obtained from a single calculation. This 
allows one to use the sum rules, at least as a theoretical tool, to understand the limitations 
with respect to the omission of core electons from the dynamics. It would be exceedingly 
challenging to ensure that the sum rules Eqs. (126|) and (127|) are obeyed in formalisms that 
require the explicit construction of the excited state spectrum. 

The violation of the sum rule Eq. f l27p in the valence particle space raises an issue that 



needs to be addressed in future work. In Ref. [3l|, it was found that the violation of 
the dipole response in TDDFT is largely justified. The dynamic contribution of the core 
electrons shifts oscillator strength down into the spectral region of valence electrons, and this 
accounts physically for the increase of the sum rule, calculated only with valence electrons 
employing the nonlocal pseudopotential in the space of valence electron excitations. Whether 
there is a related mechanism to the decrease in the I2 sum rule remains to be seen. Also, 
the pseudopotential should in principle be corrected for the gauge field associated with the 
magnetism, but that was not done here. It should be mentioned that these questions will 
also arise on calculations using the Projected Augmented Wave (PAW) method 32j, since 
this also makes the Kohn-Sham operator nonlocal. 

It was also a surprise to us to find that the MCD response may have an abnormal sign. 
This goes against the picture of an electron being excited to a higher band of orbitals and 
there undergoing circular motion in the sense given by the external magnetic field. It might 
be that strong screening destroys the simple connection to the expected classical oscillation 
picture. This raises another question for future work, to investigate in a general way the 
effects of screening on the MCD. 

Finally, we have not discussed here the sensitivity to specific density functionals. Al- 
;hough not reported, we have also carried out the Ceo calculations with the LB94 functional 



331]. This gave very similar results except for Rydberg transitions, which are consider- 



ably shifted in energy, depending on the functional. Since the observables in MCD depend 
on currents, it might also be interesting to investigate the generalized TDDFT including 
current- current interactions. 
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Appendix 

In this Appendix we apply the real-time theory to a simple Hamiltonian, a spinless 
electron in an anisotropic harmonic oscillator potential. The model is completely solvable 
making it useful in checking the coding and formulas for the TDDFT in a magnetic field. 

The Hamiltonian Hq in Eq. (1) is taken as 

^.-li + it-X (37) 

We label the eigenstates of Hq by the number of excitation quanta along each coordinate 
axis, Ifixnyriz), and we set m = e = ^ = 1 in the equations below. We take the oscillator fre- 
quencies uua to be nondegenerate, so the MCD response will only have B-type contributions. 
We first need the eigenstates in the presence of the magnetic field, expanded to first order 
in the field strength. Taking the magnetic field in the 2;-direction, the relevant perturbed 
orbitals are 

(38) 



|000,5,) = 


|000) - 


-tso\110) 


|100,S,) = 


|100)- 


-isi|010) 


|010,S,) = 


1010) - 


-isi|100) 









where 

'''' 2 (ujy±uj,)(u,ujyy/'- 
The perturbed energies of the orbitals are not needed because that perturbation is second 

order in B^. 

To get the real-time response >S'^^\ wc multiply the ground state wave function by the 
field e'^'^'^y and expand over the eigenstates, to first order in k. The required matrix elements 
of the dipole operator between ground and excited states are 

(100, BMOOO, B,) = (2a;,)-V2 (40) 
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(010,S,|/ijOOO,5,) = (2, 



-1/2 



(010,5,|/i,|000,5,) = -i^iBB,{2ujyY/^l{ujl - 



(100,5,|/ijOOO,5,) = -z/iB5,(2u;,)i/V(u;2 _ ^2) 
The initial perturbed wave function is 

\'^ky{t = 0)) = |000, B,) + ^— A^IOIO, B,) + ^sB,k-^^p^\100, B,). (41) 

The time dependence is put in by multiphng the excited states by e^**^"*. The expectation 
value of iix may then be evaluated as a function of time. The result after some simplification 
is 

^xj = -^-z — ^{cosuj - COS Uyt). (42) 

^x ^y 

The short-time response given by Eq. (27) may be verified by making a power series expan- 
sion of the cosine functions in Eq. ( 142|) . Finally, the evaluation of -Rmcd by Eq. (20) may 
be verified by carrying out the cosine Fourier transform, | /q°° dt cosut coswq^ = ^{^ ~ ^o)- 
Putting in all three magnetic moment directions, the result is 

RucB = E ^(E - ^o)-^ (43) 

It may be easily verified that -Rmcd satisfies the two sum rules, Eq. (25) and (26). 
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